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, Primordial stars are formed from a chemically pristine gas consisting of hydrogen and helium. They 
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are believed to have been born at some early epoch in the history of the Universe and to have enriched the 
interstellar medium with synthesized heavy elements before the emergence of ordinary stellar populations. 
, We study the formation of the first generation of stars in the standard cold dark matter model. We follow 

^ ' the gravitational collapse and thermal evolution of primordial gas clouds within early cosmic structures 

[ using very high-resolution, cosmological hydrodynamic simulations. Our simulation achieves a dynamic 

. range of ^ 10^" in length scale. With accurate treatment of atomic and molecular physics, it allows us to 

study the chemo-thermal evolution of primordial gas clouds to densities up to p ~ 2 x IO~*g cm~'^(nH ~ 
10 cm ) without assuming any a priori equation of state; a six orders of magnitudes improvement over 
previous three-dimensional calculations. We implement an extensive chemistry network for hydrogen, 
helium and deuterium. All the relevant atomic and molecular cooling and heating processes, including 
cooling by collision-induced continuum emission, are implemented. For calculating optically thick H2 
■ cooling at high densities, we use the Sobolev method (Sobolev 1960) and evaluate the molecular line 

Qh| opacities for a few hundred lines. We validate the accuracy of the method by performing a spherical 

Q . collapse test and comparing the results with those of accurate one-dimensional calculations that treat 

$H ' the line radiative transfer problem in a fully self-consistent manner. 

I We then perform a cosmological simulation adopting the standard ACDM model. Dense gas clumps 

d • are formed at the centers of low mass (~ 10^~^A/q) dark matter halos at redshifts z ~ 20, and they 

collapse gravitationally when the cloud mass exceeds a few hundred solar masses. To examine possible 
gas fragmentation owing to thermal instability, we compute explicitly the growth rate of isobaric per- 
turbations. We show that the cloud core does not fragment in either the low-density (nn ~ IO^°cm~^) 
I or high-density lO^^cm"'^) regimes, where gas cooling rate is increased owing to three-body molecule 

formation and collision-induced emission, respectively. The cloud core becomes marginally unstable 
against chemo-thermal instability in the low-density regime. However, since the core is already com- 
pact at that point and correspondingly the sound-crossing time as well as the free-fall time are short, 
or comparable to the perturbation growth timescale, it does not fragment. Run-away cooling simply 
leads to fast condensation of the core to form a single proto-stellar seed. We also show that the core 
remains stable against gravitational deformation and fragmentation throughout the evolution. We trace 
in Lagrangian space the gas elements that end up at the center of the cloud, and study the evolution 
of the specific angular momentum. We show that, during the final dynamical collapse, small angular 
momentum material collapses faster than the rest of the gas and selectively sinks inwards. Consequently, 
the central regions have little specific angular momentum, and rotation does not halt collapse. With 
the large physical dynamic range of our simulation, we, for the first time, obtain an accurate gas mass 
accretion rate within a IOMq innermost region around the protostar. The protostar is accreting the 
surrounding hot (T ~ 2000K) gas at a rate of M > 10~^ — lO^'^Afo/yr. From these findings we conclude 
that primordial stars formed in early cosmological halos are massive. We carry out proto-stellar evolution 
calculations using the obtained accretion rate. The resulting mass of the first star when it reaches the 
zero-age main sequece is Mzams ~ IOOM0, and less 6OM0) for substantially reduced accretion rates. 

Subject headings: cosmology: theory - early universe - stars:formation - galaxies:formation 
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1. INTRODUCTION 

The current standard theory of cosmic structure formation posits that the present-day clumpy appearance of the Uni- 
verse developed through gravitational amplification of matter density fluctuations generated in its very early history. The 
energy content of the Universe and the basic statistics of the initial density field have been determined with great accuracy 
from recent observations of the cosmic microwave background (Spcrgel ct al. 2006), large-scale structure (Tegmark ot al. 
2004; Cole et al. 2005), and distant supernovae (Riess et al. 2004; Astier et al. 2006). It has become possible to make 
accurate predictions for the formation and nonlinear growth of large-scale structure within this framework. 

The standard model based on cold dark matter (CDM) predicts that the mass variance after the epoch of matter- 
radiation equality has progressively larger amplitudes on smaller mass scales, indicating that objects form hierarchically, 
with smaller ones originating first. Although the fluctuation amplitudes on scales much smaller than galactic scales are 
not known directly from observations, most of the viable candidates for the dark matter predict nearly scale-invariant 
fluctuations down to stellar or even planetary masses. Such CDM models generically argue that the first cosmological 
objects are formed in low mass (~ IO^Mq) dark halos at high redshifts (z ^ 20) when primordial gas condenses via 
molecular hydrogen cooling (Couchman & Rees 1986; Haiman, Thoul & Loeb 1996; Tegmark et al. 1997; Abel, Bryan & 
Norman 2002; Yoshida et al. 2003; see Bromm & Larson 2004 for a review). Hierarchical structure formation eventually 
leads to the emergence of a population of early generation stars/galaxies which terminate the Cosmic Dark Ages by 
emitting the first light (Madau 2000; Cen 2003; Miralda-Escude 2003; Sokasian et al. 2003, 2004; Reed et al. 2005). 

The first stars are also believed to be the first source of heavy element production (those heavier than lithium). Heavy 
elements must have been processed in massive stars and expelled into the interstellar medium by supernovae, to enable 
the formation of ordinary stellar populations. Early metal-enrichment of the low density intergalactic medium is also 
suggested by observations of the Lyman-a forest (Songaila 2001, 2005; Schaye et al. 2003). Interestingly, the recent 
discovery of extremely metal-poor stars (Christlieb et al. 2002; Frebel et al. 2005) suggests that relics from very early 
generation stars exist even today in our own Galaxy. Stars with such low metallicities provide valuable information on 
the first stars that likely polluted the interstellar medium (Iwamoto ct al. 2005; Tumlinson 2006). 

The study of primordial star formation has a long history. The formation of the first cosmological objects via gas 
condensation by molecular hydrogen cooling has been studied for many years (e.g. Saslaw & Zipoy 1967; Peebles & 
Dicke 1968; Matsuda, Sato, & Takeda 1969; Kashlinsky & Rees 1983). In the context of star formation, the evolution 
of a collapsing primordial gas cloud has also been studied extensively (Matsuda, Sato, & Takeda 1969; Yoneyama 1972; 
Carlberg 1981). Palla, Salpeter & Stabler (1983) used a one-zone model to follow gas collapse to very high densities. 
They identified a minimum Jeans mass scale of 0.1 solar mass from the thermal evolution of the cloud core, suggesting 
the formation of low-mass primordial stars. One-dimensional hydrodynamic simulations of spherical gas collapse were 
performed by a number of researchers (Villere & Bodenheimer 1987; Omukai & Nishi 1998; Ripamonti et al. 2002). 
Omukai & Nishi (1998) included a detailed treatment of all the relevant chemistry and radiative processes and thus were 
able to provide accurate results on the thermal evolution of a collapsing primordial gas cloud up to stellar densities. These 
authors found that, while the evolution of a spherical primordial gas cloud proceeds in a roughly self-similar manner, there 
are a number of differences in the thermal evolution from that of present-day, metal- and dust-enriched gas clouds. 

Recently, three-dimensional hydrodynamic calculations were performed by several groups. Abel et al. (2000; 2002) and 
Bromm et al. (1999; 2002) studied the formation and fragmentation of primordial gas clouds in CDM models. From 
various analyses of the simulation results, these authors conclude that primordial stars (often called "Population III") 
that form in dark matter matter halos are rather massive. Yoshida et al. (2003) examined the statistical properties of 
primordial star-forming clouds from a large sample of early mini-halos located in their cosmological simulations. The 
chemo-thermal evohition of the primordial gas clouds is rather complex, being coupled with the dynamical evolution of 
the assembly of dark matter halos (Yoshida et al. 2003; Gao et al. 2005), and simple analytic models based on crude 
assumptions almost always fail in predicting the evolution of primordial gas in CDM halos. 

The three-dimensional calculations of Abel et al. (2002; hereafter, ABN02) were able to follow the gas evolution until 
a small (~ IMq) fully molecular core formed. However, since they did not include the effect of the opacity of the 
cloud core, the radiative cooling rate is over-estimated at the cloud center. Thus their results for late stage evolution, 
including the important problem of fragmentation and mass accretion, are uncertain. Bromm & Loeb (2004) made similar 
approximations and also continued their simulation by creating a sink particle to study the late time accretion phase. 
None of these simulations reproduce the correct density, temperature, and velocity profiles in the vicinity the protostar, 
which are among the most important quantities for proto-stellar evolution. A critical technique we describe in the present 
paper is compiitation of molecular line opacities. With the implementation of optically thick line cooling, the gas evolution 
can be followed to much higher densities than in the previous studies. In principle, evaluating the net line cooling rate 
at such high densities involves costly line radiative transfer. We circumvent this by using a local Sobolev method which 
employs local velocity information as well as knowledge of the density and temperature. We show that the method works 
well in problems of collapsing gas clouds, in terms of computation of radiative cooling rates and resulting density and 
temperature structiire. Wc apply this technique to cosmological simulations. 

An important, outstanding issue is to determine whether or not a gas cloud fragments into multiple objects during its 
collapse. If it does fragment, it could yield multiple stars, suggesting the possibility of star-cluster formation including 
low-mass stars (Sabano & Yoshii 1977; Silk 1983). To the contrary, if the gas cloud remains stable against fragmentation 
throughout its evolution, it will lead to the formation of a single proto-stellar seed, surrounded by a large amount of 
infalling gas - a typical condition for a proposed scenario of the formation of massive stars (Larson & Starrfleld 1971; 
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sec the review by Larson 2003). We show that the gas cloud formed at the center of a high-redshift halo remains stable 
against fragmentation, to form a single proto-stellar seed. We also find that the rate of gas accretion is very high. The 
mass accretion rate onto a protostar largely affects the final stellar mass (Omukai & Palla 2001; 2003). Stahler, Palla, & 
Salpeter (1986a, b) suggested that large accretion rates (> 1O~^M0 yr~^) are expected for primordial protostars because 
of high gas temperatures in the infalling surrounding gas. For a better evaluation of the mass accretion rate, Omukai & 
Nishi (1998) use a self-similar solution to derive a time-dependent rate, while Tan & McKee (2004) consider effects of 
rotation in accretion dynamics around proto-stellar disks. We address this issue using fully self-consistent cosmological 
simulations. Finally, we carry out proto-stellar evolution calculations and determine the resulting stellar mass on the 
assumption that the true gas mass accretion rate and its time-dependence is close to the instantaneous accretion rate at 
the protostar formation epoch. 

The rest of the paper is organized as follows. In Section 2, we describe the chemistry network used for our simulations. 
The relevant cooling and heating processes are described in Section 3. There, we also provide details of numerical 
implementation. In Section 4, we present the results of a test problem of a spherically collapsing gas. We compare our 
results to previous work using one-dimensional codes. We present the results of our ACDM cosmological simulations in 
Section 5. Implications of the main conclusions are discussed in Section 6. 

2. PRIMORDIAL GAS CHEMISTRY 

The reaction network for hydrogen and helium (c~, H, H+, He, He+, He^"'", II2, Hj, H~) is largely based on that of 
Abel et al. (1997) and Galli & Palla (1998; hereafter, GP98). The rate coefficients are summarized in Table 1. In this 
section, we describe updates in detail. 

2.1. Molecule formation 

In a low density primordial gas, the main formation path for hydrogen molecules is via the H~ channel: 



H + e^B.'+hi' (reaction 7), (1) 

H -I- H" ^ H2 -I- e (reaction 8), (2) 

where electrons are used as catalysts. H2 formation via the 112" channel (reaction 9, 10) is dominant at very high redshifts 
z > 200 where cosmic microwave background photons destroy H~ by photo-detachment (sec also a recent calculation of 
Hirata & Padmanabhan 2006, which shows the importance of reactions involving HeH-l- ions in H2 production at such 
high redshifts.) There is considerable variation among published experimental data and theoretical calculations for the 
reaction rate /cg. However, uncertainty in this rate will have little effect on the H2 abundance, because the reaction is 
always much faster than the formation of H~. We use the fit by GP98 to the calculations of Launay et al. (1991) for kg- 
The rate coefficients for mutual neutralization 



H--FH+^H2 (reaction 18), (3) 

are also uncertain, particularly at low temperatures (Glover et al. 2006). This reaction affects the H2 abundance 

significantly only when the ionization fraction is greater than ^ 0.01 (Palla & Zinnecker 1987). In the calculations 
presented below, the ionization fraction is always smaller than 10^'^, and hence the uncertainty in the reaction rate is 
unimportant. We use the fit given in GP98. We note that a better determination of this rate is clearly needed for 
calculations of, for example, the evolution of cooling gas in relic Hii regions (Oh & Haiman 2002; Glover et al. 2006; 
Yoshida 2006; Susa & Umemura 2006). 

2.2. Three-body reactions 
At high densities {nu > lO^cm"'^), hydrogen molecules are formed by the following rapid reactions 

H + H + H->H2+H (reaction 22), (4) 

and 

H -h H H2 ^ 2H2 (reaction 23). (5) 
We use the reaction coefficients as given in Palla, Salpeter, & Stahler (1983): 

A;22 = 5.5 X 10-29 cm^s'S (6) 

and 

k23 = 6.875 X 10"^° T"^ cm^s'^ (7) 

Note that the net formation rate scales with the cube of density. The three-body reactions are very efficient in converting 
hydrogen atoms to molecules at densities n > 10*cm~^. 

Krstic, Janev & Shultz (2003) recently revised reaction rates for three-body reactions involving protons, 

H + H + H+^H2 + H+, (8) 
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H + H + H+ ^ H+ + H, (9) 

and found the rates are as large as those for reactfons (22), (23). Since the ionization fraction is extremely small in the 
primordial gas at high densities and low temperatures (T < 3000 K), which is the range of our interest, these reactions are 
unimportant in the calculations presented in what follows. We note, however, that when the ionization fraction (hence 
proton abundance) is large in the proto-stellar gas cloud in the final adiabatic phase, or under strong ultra-violet radiation, 
inclusion of the reactions could be important. 

2.3. Collisional dissociation of H2 
In a neutral primordial gas, H2 molecules are destroyed mostly by collisions with H atoms and H2 molecules: 

Ha + H^SH (reaction 11), (10) 

H2 + H2 ^ H + H + H2 (reaction 25) (11) 

We use the collision cross-sections of Martin et al. (1996) for reaction 11 that include dissociation from excited levels at 
high densities. For dissociation by collisions with H2 (reaction 25), we use the high density limit of Palla et al. (1983) 



fc25 = 8.125 X 10-« T-V2exp ( _52000\ 



, 6000\ 
J X ^1.0 - exp —j 



(12) 



T 

Dissociation by collisions with helium atoms 

H2 + He ^ H + H + He (13) 

and its reverse reaction (three-body reaction involving helium) have much smaller reaction rates than those involving 
collisions with hydrogen (Dove et al. 1987), and thus we ignore these processes. 

We include the charge exchange reaction between H and He (reaction 14, 15, see Glover & Brand 2003) and update 
the reaction rate for H + HJ H+ + H2 (reaction 12) with the new fit given by Savin et al. (2004a, b). These processes 
become important in calculations of, for example, the evolution of relic Hii regions (Yoshida 2006), but do not affect the 
calculations in the present paper. 



2.4. The adiabatic index 
We write the equation of state for a gas consisting of atoms and molecules as 

kBT 



P = 



limp 



-P = (7 - 



(14) 



where the adiabatic index is computed from 




(15) 



The summation is taken over all chemical species including helium. For the adiabatic exponent for H2 , we use the formula 
in Landau & Lifshitz (1984), 



7H2 



1 



5 + 2a:^ 



(e^ - 1)2 



(16) 



where the last term with x = 6100K/r accounts for vibrational degrees of freedom of hydrogen molecules. At very high 
pressures (P » 10*dyne cm~^), various non-ideal gas effects become important and the equation of state needs to be 
modified (Saumon, Chabrier, & van Horn 1995; Ripamonti et al. 2002), but this is outside the range of gas phases we 
consider here. 



3. COOLING AND HEATING RATES 

Radiative cooling processes owing to excitation, ionization, and recombination of hydrogen and helium atoms are well 
established. We refer the readers to Cen (1992)^ , Fukugita & Kawasaki (1994), and Abel et al. (1997). We note that these 
cooling processes are generally unimportant in studies of the evolution of a neutral primordial gas. Below we describe 
more relevant molecular cooling and chemical heating processes in detail. 

^The recombination rate for He++ in Cen (1992) has an incorrect temperature dependence. See Fukugita Sz KawasaJji (1994) for the correct 
scaling. 



5 



3.1. Cooling by H2 and HD molecules 
We use the cooling rate of Galli & Palla (1998) for H2 line cooling in the low density limit: 

Kn^n ^ 0) = dex[-103.0 + 97.59 logT - 48.05(logT)2 + 10.80(logr)3 - 0.9032(logT)'*] 



(17) 



At high densities, formation of H2 molecules takes place and molecular fraction increases rapidly, making H2 line cooling a 
dominant cooling process. At high densities, all energy levels are populated according to local thermodynamic equilibrium 
(LTE). An accurate fit for the optically thin cooling rate is given by HoUenbach & McKee (1979): 



A 



LTE 



where 



and 



_ 9.5 X 10-22 r3-76 
Arnt — 



1 + 0.12 r|i 



exp 



= A,c 



0.13 



A. 



vibi 



+ 3 X 10-24 cxp 



A. 



6.7 X 10-^^ exp 



5.86 



1.6 X 10"^** exp 



11.7 

~t7 



0.51 



T3 = 



T 



lOOOK 



The low density hmit and the high density limit are bridged as 



Ah. = 



Ah. (LTE) 

:r/nH ' 



1 



where 



^ Ah, (LTE) 



(18) 
(19) 

(20) 

(21) 
(22) 



For HD cooling, we use the formulation of Flower et al. (2000) ^ . We show the cooling rates in Fig. There, we assume 
a fractional abundance of [HD/H2] — . Cooling by HD molecules is important only at low temperatures {T < 200K) 
and low densities (tt-h < 10^); otherwise H2 cooling dominates. We have carried out a spherical collapse simulation (see 
Section^ with full deuterium chemistry for five species (D, D+, D~, HD, HD+), using the extensive reaction network of 
Nakamura & Umemura (2002). We find that the fractional ratio of [HD/H2] has a maximum of ^ 10~^ at nn ^ lO^ cm^^, 
where the gas temperature is about 200 K (see Fig. Oand also Omukai et al. 2005). Therefore, HD cooling is unimportant 
for the thermal evolution of initially neutral primordial gas clouds. The role of HD cooling in partially ionized gas will be 
studied elsewhere (Yoshida et al. 2006, in preparation). 



10-^ 



o 




10 100 1000 10000 

T[K] 

Fig. 1. — Radiative cooling rates per molecule as a function of temperature. We show the low-density limit for H2 (solid line) and HD 
(dot-dashed line) assuming njj = 1 cm"'^ and the relative abundance HD/H2 = 0.001. The dashed line is the high-density limit of the H2 
cooling function. 

^The HD cooling function of Flower et al. does not include transitions between high vibrational levels but the contribution to the total 
cooling rate at the relevant densities are unimportant, as shown by Lipovka et al. (2005). 
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3.2. Chemical cooling and heating 
We include the heat gain associated with the formation of hydrogen molecules via three-body reactions as 

d?iH2 



rH2,3b — 



dt 



(23) 



where — 4.48eV is the molecular binding energy. We also account for heat loss owing to coUisional dissociation (Section 
I2.3|l in the same manner. The two-body formation process of H2 molecules, either via (reaction 8) or (reaction 
10), is exothermic, but deposits its energy into H2 through rotational and vibrational excitation. These reactions can 
be important sources of heat at high densities where coUisional de-excitation removes the excitation energy. Following 
HoUenbach & McKee (1979) and Shapiro & Kang (1987), we calculate the heating rates as 



r(H") fc8nH-nH[3.53(l + ricr/nn)"^] eVcm-3s-\ 
r(H+) = fcionH+7iH[1.83(l + n„/nur^] cVcm^^s^S 
where the critical density is given by 



= lOT 



-1/2 



1.6yHexp[-(400/T)2] + 1.4 yn. exp - 



T 



12000 \ 
1200 j 



cm 



(24) 
(25) 

(26) 



and 2/H and 2/H2 are the number fractions of hydrogen atoms and hydrogen molecules, respectively. 

3.3. Optically thick H2 line cooling 

When the gas density and the molecular fraction are high, the cloud becomes opaque to molecular lines and then H2 
line cooling becomes inefficient. The net cooling rate can be expressed as 

Ah 2, thick — 

u.l 



(27) 



where n„ is the population density of hydrogen molecules in the upper energy level u, A^i is the Einstein coefficient for 
spontaneous transition, /3csc,ui is the probability for an emitted line photon to escape without absorption, and hvui — ^E^i 
is the energy difference between the two levels. In the relevant temperature range T ~ 1000 — 2000K, we consider rotational 
levels from J = to 20, and vibrational levels v = 0,1, 2. We use the radiative transition rates Aui for these levels from 
Turner, Kirby-Docken, Dalgarno (1977). We calculate the vibrational energies E{v,J) following Borysow, Frommhold 
& Moraldi (1989). Note that, although the net cooling rate can be directly obtained from equation H27|l . we use it to 
evaluate the reduction factor /red = Athick/Athin- For computational efficiency, we compute the cooling rate as a function 
of temperature using equations (|18|I - H20|I . and multiply the rate by /red when nn > lO^cm"'^. 

In order to calculate the escape probability, we first evaluate the opacity for each molecular line as follows. The 
absorption coefficient for a transition from /-level to w-level is 



in 



niBi, 



1 — exp 



kT 



where ni is the number density of hydrogen molecules in level I, Biu is Einstein's B-coefl[icient, and 
We approximate the line profile by a Gaussian function 



(28) 



{v) is the line profile. 



1 



exp 



(29) 



with the thermal Doppler width /S.vd — {vq / c) / kT / mY{. We then compute the opacity at the line center as 

Tin = aiuL, (30) 

where L is the characteristic length scale, which is approximately the cloud core size. Since the absorption coefficients are 
computed in a straightforward manner, although somewhat costly, the remaining key task is the evaluation of the length 
scale L. 

It is important to recall that calculation of the photon escape probability is essentially a line transfer problem, for 
which, in principle, demanding radiative transfer calculations are needed to obtain accurate results. However, by noting 
that the important quantity we need is the effective gas cooling rate, we can formulate a reasonable and well-motivated 
approximation. To this end, we decided to use the escape probability method that is widely used in the study of stellar 
winds and planetary nebulae (Castor 1970; Goldreich & Kwan 1974). Consider a photon traveling in a gas in one direction 
along which a constant velocity gradient exists. The escape of photons is greatly enhanced by the presence of macroscopic 
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velocity fields, i.e. large velocity gradients; when an emitted photon has traveled one Sobolev length to the point where 
the profile is Doppler-shifted by one characteristic width (line width), it can travel unimpeded and will escape from the 
local neighborhood. 

We calculate the Sobolev length along a line-of-sight as 

^ ^thermal /'Ql\ 

" \dVr/dr\ ' ^^^^ 



where Wthormai = \/ kT/niii is the thermal velocity of H2 molecules, and Vr is the fluid velocity in the direction. A suitable 
angle-average must be computed in order to obtain the net escape probability. For a spherical cloud, the escape probability 
is given by 

1 — cxp(— r) 

/3esc = , (32) 

T 

with T — aLr (Castor 1970; de Jong, Dalgarno & Chu 1975). In three-dimensional calculations, in which we do not 
know the exact geometry and alignment of an object, we compute the Sobolev length and the escape probability in three 
arbitrary orthogonal directions and take the mean as 

p=^l±A±Jk, (33) 
3 

We first test the implementation in this section, and the overall accuracy of the method is examined in Section ^ 
Consider a homologous contracting sphere with size R — 0.01 pc. The velocity at radius r is given by 

v{r)=rV/R, (34) 

where = 10 km/sec is the radial (inward) velocity at R. We set the gas temperature to be a constant at T = IGOGK 
and assume the gas is fully molecular with nH2 — 5.0 x 10^°cm~^. In this case, the optical depth for an emitted photon 
(owing to a level transition it — > /) is a constant 



Tlu — —. nil^iu 



1 — exp 



kT 



and the escape probability is 



1 - exp(-T;„) 

Pesc,i« = • (3d) 

Tlu 

We run a three-dimensional simulation with the above set-up using a half million gas particles and compute the escape 
probability (equation \3'S\ ) at random points in the simulation region. For illustrative purpose, we consider only the 
rotational transition J — 6 4(i> = 0), which is one of the strongest lines at T '--^ lOOOK. Putting all the numerical 
values into equation (|35|l . we find t^q — 3.61 and P^q = 0.269. From the test simulation, we obtained a mean value of 
P = 0.2692 within r < 0.75i? and a maximum deviation of less than 3 %. We did not include the outer parts of the sphere 
(r > 0.75i?), where the density estimate is affected by the boundary. These results are quite satisfactory, and assure us 
that the numerical implementation is done correctly. We discuss the overall accuracy of this method in Section 0] where 
we perform a gas collapse simulation with density and temperature gradients. We emphasize that simple local estimates, 
for example, using only local densities and temperatures, do not provide correct values for the optically-thick line cooling 
rate, because the escape probability of line photons is greatly affected by local velocity gradients which are highly variable 
in both space and time. 

The Sobolev method can also be extended to calculation of the self-shielding factor against photo-dissociation by 
Lyman- Werner photons. The feedback effects of far ultra-violet radiation is often modeled either by assuming the gas is 
optically thin to the photo-dissociating photons (Machacek, Bryan & Abel 2001; Ricotti et al. 2002), or by employing 
the formula of Draine & Bertoldi (1996) in cases with simple geometry (Omukai 2001). While accurate results can be 
obtained by solving radiative transfer for a number of lines, as done by, e.g.. Glover & Brand (2001) in one-dimension, 
it will be useful to devise an approximation. The first such attempts of accounting for self-shielding have been already 
made, although in a much simplified manner, by Yoshida et al. (2003) and Glover et al. (2006). Our method (or its 
simple extension) may provide better solutions than these previous works. 

3.4. Cooling by collision-induced emission 

At densities greater than nn ^ lO-'^^cm^'^, hydrogen molecules collide frequently and collision pairs can generate an 
induced electric dipole, through which either molecule make an energy transition by emitting a photon. This process is 
known as collision-induced emission (CIE), the opposite process to collision- induced absorption: 

H2(v, J)+H2 ^H2(w', J')+H2 + /ii^, (37) 
H2(w, J)-hHe^H2(w', J') + He-H/ii/. (38) 
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Fig. 2. — C ooling rate per molecule owing to collision-induced emission as a function of temperature. The solid points are calculated from 
equation We consider only the dominant H2-H2 and H2-He collisions assuming hydrogen is all in molecules. The solid line is our fit 

equation 1401 . We assume nH2 = lO^^cm"'' in this plot. 

This process yields very complex spectra, and have an essentially continuum appearance. For H2 -H2 collisions, only 
fundamental vibrational bands can be noticed as smooth bumps in the spectra (see Formmhold 1994 for full account of 
this process). Thus we need to calculate the total emissivity by integrating the contribution from each transition: 



VCIE = 



^2 crciE n(H2)exp ( -^^ 



hi> 



(39) 



We use the cross-sections of Jorgensen et al. (2000) and Borysow, Jorgensen, & Fu (2001) and Borysow (2002). In 
practice, for simulations of a collapsing primordial gas, CIE cooling is important in the phase where the gas is nearly fully 
molecular (Omukai 2001). We consider H2 -H2 and H2 -He collisions, where the latter gives a minor contribution to the 
cooling rate. In practice, we use a simple fit for the total cooling rate as a function of temperature as 



AciE =dex[-116.6 + 96.34 X logT 

- 47.153 X (logT)2 -f- 10.744 x (logr)^ - 0.916 x (logT)"*] 



(40) 



In Fig. 13 we compare the fit with the accurate cooling rate obtained directly from the cross-sections. Clearly, equation 
(|40|l provides an excellent fit over the plotted temperature range. 

3.5. Limitation of the simulations 

We successfully implemented the relevant atomic and molecular physics to follow the gas evolution to very high densities, 
riH ~ lO^^cm"^. What we have included is sufficient for the present work, but in principle we will need to implement 
additional physics to probe higher densities. We here briefly describe necessary modifications in order to follow the 
evolution of the proto-stellar gas further, up to stellar densities. First, for densities much larger than lO^^cm"^, the 
reaction time scale becomes significantly shorter than the dynamical time. Then, we can use equilibrium chemistry 
which would actually simplify our handling of chemistry evolution. The species abundances can be determined from the 
Saha-Boltzmann equations, for example. Second, for rifj > lO^^cm"^, the gas cloud becomes opaque even to continuum 
emission, and then we need to evaluate continuum opacity. We have already implemented the calculation of the continuum 
opacity using the table of the Planck opacity for primordial gases (Lenzuni, Chernoff, & Salpeter 1991; Mayer & Duschl 
2005). Third, when dissociation of hydrogen molecules is completed at T ~ 5000 K, there will be no further mechanisms 
that enable the gas to lose its thermal energy, and then the gas temperature increases following the so-called adiabatic 
track. The equation of state in the high pressure regime must be modified to account for non-ideal gas effects (Saumon 
et al. 1995; Ripamonti et al. 2002). In future work, we pursue simulating the formation of protostars by employing much 
higher mass resolution and all these necessary physics. 

3.6. Time integration 

We incorporate the chemistry solver in the parallel A^-body/Smoothcd Particle Hydrodynamics (SPH) code GADGET- 
2 (Springel 2005) in the following manner. The code computes gravitational and hydrodynamic forces and updates the 
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particles positions and velocities. All the force terms as well as other quantities are evaluated in double precision in our 
simulations. For gas particles, the code also updates the density and specific entropy. In particular, we employ the fully 
conservative formulation of SPH (Springel & Hernquist 2002), which maintains energy and entropy conservation even 
when smoothing lengths vary adaptively (Hernquist 1993). This refinement to the SPH method is crucial for accurately 
describing the evolution of gas over the large dynamic range required to follow primordial star formation in a cosmological 
context. At the end of one time step, we evolve the fractions of the chemical species, and add an extra heating (cooling) 
term which arises from chemical reactions. 

We use a backward difference formula (Anninos et al. 1997) 

For each gas particle, we supplement the usual Courant condition with two additional constraints so that the time step 
does not exceed the characteristic cooling time and the characteristic chemical reaction time. We monitor the cooling 
time 

T 

Atcooi = etoi^ (42) 

and use the rate of change of the electron number density and that of neutral hydrogen to measure the characteristic 
chemical reaction time 

Atchcm = etoi min{^ , }. (43) 

Setting etoi = 0.1 suffices for following the evolution of low density gas {n < lO^cm"^). At high densities, particularly 
at n ^ 10^ — 10^^cm~'^, rapid reactions release a significant amount of heat. After some experiments, we found that the 
chemistry part becomes numerically unstable with etoi = 0.1 and concluded that setting etoi = 0.01 allows stable time 
integration at n > lO^cm"'^. 

4. SPHERICAL COLLAPSE TEST 

In this section, we test our numerical techniques using a spherical collapse problem. Wc follow the evolution of primordial 
gas in a dark matter halo by setting up a gas sphere embedded in a NFW (Navarro, Frenk, White 1997) potential 

P'^^^ = / / N/1 ^% / ^2^ ' (44) 

(r/rs)(l + [r/rsy) 

where rs,Ps are scale radius and density, respectively. The initial gas density is set to be an isothermal f3 profile 

= [i + (,%V]3,/2 - (45) 

We are interested in the gas evolution after the gas cloud becomes self-gravitating, and so details of the initial density 
profile do not matter. For simplicity, we set /3 = 1. The halo mass is set to be 5 x 10^ Mq and we assume the baryon 
fraction to be 0.05. We distribute 4 million particles according to equation H45|l and evolve the system. 

Fig. 131 shows the distribution of gas in a temperature-density phase plane when the central density is 5 x lO^^cm"^. 
In the figure, characteristic features are marked as regions A-G. The bottom panel in Fig. O shows the corresponding 
molecular fraction distribution. All the features are explained as closely related to the thermal evolution. See the caption 
for a brief explanation. The overall evolution of the central gas cloud after it undergoes a run-away collapse is consistent 
with the spherically symmetric calculation of Omukai & Nishi (1998; hereafter, ON98). We have checked the radial profiles 
of density, temperature, velocity and molecular fraction. These quantities are quite similar to the late time evolution of 
the ON98 calculation until the central gas density reaches ~ lO^^cm"^ (up to the fourth output in Fig. 1 of ON98). 
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Fig. 3. — (Top) Gas distribution in the temperature-density phase space in a spherical collapse problem. The indicated characteristic 
features are explained as follows; (A) gas temperature reaches > lOOOK by virialization, and hydrogen molecules are formed by two-body 

processes, (B) molecular hydrogen cooling brings the gas temperature down to 200 K, (C) the H2 cooling rate saturates and becomes close 
to the density-independent, LTE value, (D) three-body reactions kick in and the gas becomes fully molecular, (E) the line cooling rate 
decreases as the density increases because of the cloud's opacity, (F) collision-induced emission becomes a dominant cooling process, and (G) 
H2 dissociation begins at T ~ 2000K. (Bottom) The molecular fraction /h2 of the gas. The increase in the fraction at A, D, a plateau at 
C— >D, and the temporal decrease owing to dissociation at G are clearly seen in this plot. 

Previous three-dimensional simulations of primordial gas cloud formation were hampered by the complexity of calculat- 
ing line opacities and the reduction of the resulting cooling rate. The maximum resolution, in terms of gas density, achieved 
in these simulations was thus ^ 10^*^ cm"'^, where the assumption of optically thin cooling breaks down. With the novel 
technique described in Section 13.31 our simulations can follow the evolution of a primordial gas cloud to rin ~ lO^^cm"^, 
nearly six orders of magnitude greater than previous three-dimensional calculations reliably probed. To study the detailed 
evolution of a proto-stellar "seed" beyond njj ''^ lO^^cm"^ , we would need to implement a few more physical processes, 
as explained in Section 3.5. 

An important quantity we measure is the optically-thick line cooling rate, which serves as a critical check of our 
numerical implementation. Fig. ^ shows the normalized H2 line cooling rate against local density. We use an output at 
the time when the central density is ric — lO^^cm^'^. In the figure, we compare our simulation results with those from 
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the full radiative transfer calculations of ON98 (open diamonds). Clearly our method works very well. The steepening 
of the slope at n > lO^^cm"^, owing to the velocity change where infalling gas settles gradually onto the center, is well- 
reproduced. We emphasize that the level of agreement shown in Fig. 0]can be achieved only if all of the local densities 
(of chemical species), temperatures, and velocities are reproduced correctly. 
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Fig. 4. — The cooling efficiency defined by / = Atjjick/Athin- The open squares are the results from the one-dimensional calculation of 
Omukai & Nishi (1998). We plot the efficiency as a function of local density at the time when the central density is ra = lO^'^cm"''. 

5. COSMOLOGICAL SIMULATIONS 

It is important to study the formation of primordial stars using a proper set-up, starting from realistic initial conditions. 
To this end, we have run a cosmological simulation adopting the concordance ACDM cosmology with matter density 
ftm ~ 0.26, baryon density fib = 0.04, cosmological constant f^A = 0.7, and expansion rate at the present time Hq ~ 70km 
s^^ Mpc~^. The power spectra of the initial density fluctuations for baryons and dark matter are calculated by the 
Boltzmann code of Sugiyama (1995). In the examples discussed below, the power spectra were normalized to erg = 0.9, 
consistent with the first year WMAP results, but slightly larger than the recent estimate of Spergel et al. (2006), based 
on their third year data-set. Everything else being equal, a smaller value of as will shift structure formation to slightly 
lower redshifts. The main purpose of our present investigation is to study the physics of gas collapse to high densities 
and fragmentation in dark matter halos. These processes will be unaffected by a shift in ag, although the cosmic time 
corresponding to the evolution we model would be different, as would the abundance of the halos and stars at a given 
redshift. 

In order to avoid spurious clumping in the initial particle set-up, we use "glass" particle distributions (White 1996). 
Further details on the initial conditions are found in Yoshida, Sugiyama & Hernquist (2003) . We use a simulation volume 
of 0.3 Mpc on a side. Starting from a low resolution simulation, we select the most massive halo at z = 15 in the volume 
and apply a hierarchical zoom-in procedure (e.g. Navarro & White 1993; Tormen et al. 1997; Gao et al. 2005) to the 
region surrounding the halo, so that progressively higher mass resolutions are realized. In the highest resolution part 
of the the initial conditions, the dark matter particle mass is 0.0944 Mq and that of gas particles is 0.0145 Mq. In 
order to achieve an even higher mass resolution, we progressively refine the gas particles using the method of Kitsionas 
& Whitworth (2002) and Bromm & Loeb (2003) as the gas cloud collapses. We do this refinement so that a local Jeans 
length is always resolved by fives times the local SPH smoothing length. With this technique, we achieve a mass resolution 
of nip — 6OM0, where denotes the Earth mass, at the last output time. For the study of cloud fragmentation at 
intermediate densities rijj = 10^ — lO^^cm"^, we carry out two simulations; one with the refinement technique and the 
other without it, to verify that temporal perturbations caused by refinement do not affect the results presented in Section 
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Fig. 5. — Projected density distribution for our cosmological simulations at a redshift z f» 19. The physical side-length is indicated in each 
panel. 
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Fig. 6. — Radial profiles for density, temperature, molecular fraction, and infall velocity at a redshift z 19. The density is plotted as a 
function of distance from the center, whereas the other three quantities are plotted as a function of enclosed gas mass. The density profile is 
close to the power-law oc R~^'^. The characteristic features in the temperature profile (marked A-G) can be explained as in Fig. |3](see also 
text). 



5.1. Thermal evolution of the proto-stellar gas cloud 

In our ACDM simulation, the first gas cloud is formed at the center of a dark matter halo with a mass of 6 x 10^ Mq. 
The halo is located at the intersection of filamentary structures at a high-density peak. Fig. [S] shows the projected gas 
density around the proto-stellar gas cloud. At the output time, the density of the cloud core has reached '--^ 3 x lO^^cm"^. 
The four panels in Fig. |Slshow progressively zoomed-in views of the central regions, with the physical dimension indicated 
in each panel (clockwise from 200 parsec to 100 AU). In the lower-left panel, we compare the size of the cloud core 
with a solid bar which indicates a length of 5 astronomical units! The central part in the top-right panel has a mass of 
~ 300Mq with an approximate diameter of ~ Ipc; this region is self-gravitating and undergoing run-away collapse. In 
the bottom- right panel, the central darkest area is nearly fully molecular with a mass of about IM©. Finally the darkest 
region in the bottom-left panel has a mass of ^ O.OIMq, which is completely opaque to molecular lines and is cooling by 
CIE. 

Fig. 1^1 shows the radial profiles of density, temperature, velocity, and molecular fraction at the final output time. 
Except for the density profile which is shown as a function of radius (note the x-axis is in astronomical units), the mass 
coordinate (enclosed gas mass) is used for the horizontal axis. Although the plotted profiles are for a single output 
time, the gas evolution can indeed be easily inferred from them, because the collapse time scales as oc l/\/nj leaving the 
profiles of the outer region almost unchanged while the central parts collapse further. The density profile evolves self- 
similarly approximately as a power-law with n cx r"^'^ (dashed line in the top-left panel), being consistent with previous 
one-dimensional and three-dimensional simulations. 

The characteristic features of the temperature profile are understood by various physical processes as we discussed 
already in the spherical collapse simulation. Virialization brings the gas temperature to above 1000 K, and hydrogen 
molecules are formed by two-body reactions (point A in the top- right panel of Fig. EJ. The cloud cools by H2 cooling 
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(B) and the gas temperature is lowered to about 200 K (C). There, the gas density is ^ lO'^cm^'^, beyond which the H2 
rotational level population approaches that of LTE. Then, the cooling rate per molecule becomes density-independent, 
and radiative cooling is much less efficient than in the low density limit. This is the so-called loitering regime (Bromm 
et al. 2002), and further collapse is induced dynamically when the cloud becomes Jeans- unstable. To see the onset of 
run-away collapse, we compare the enclosed gas mass with the locally estimated Bonnor-Ebert mass (Bonnor 1956; Ebert 
1955): 

M - '^^''^ 



G3/2pV2 ' 

« 20AfQr^/2n~i/V"^7^ , (46) 

where mi is the first maximum mass of the solution for the isothermal Lane-Emden equation (see, e.g. Stabler & Palla 
2004), and /z and 7 denote the mean molecular weight and adiabatic index, respectively. We approximate the external 
pressure by its local value taken from the radial density and temperature profiles. (Note that Mbe evaluated in this 
manner is essentially the same as the local Jeans mass.) We find that the enclosed gas mass exceeds the Bonnor-Ebert 
mass at M ~ 200Mq. This is the characteristic mass of the collapsing gas cloud. We see in Fig. Elthat the temperature 
rises below this mass scale, showing clearly the onset of collapse and heating by contraction (Point C). The 2OOM0 cloud 
is now self-gravitating, being essentially decoupled dynamically from the host dark matter halo. 

Three-body reactions convert most of hydrogen into molecules (D in the top- right panel). As the molecular fraction 
profile (bottom-left panel) shows, the central 1-2M0 becomes fully molecular through these rapid reactions. Cooling by 
H2 lines soon saturates because of gas opacity (E). When the core density reaches ~ lO^^cm^'^, H2 line cooling becomes 
very inefficient, but CIE cooling kicks in (region F). At higher densities, the gas temperature is maintained at ~ 2000K, 
whereby dissociation of hydrogen molecules effectively absorbs heat input from contraction and dissipation of gravitational 
energy (G). Dissociation already commenced in the innermost region at the plotted output time, as is clearly seen in the 
molecular fraction profile at Mcnc < O.IMq. 

5.2. Thermal instability and fragmentation 

It is an outstanding question whether or not primordial gas clouds are subject to thermal instability and fragment into 
multiple objects. Fragmentation could be manifested in various forms because it depends on the cloud geometry, thermal 
history, interaction with ambient matter, etc. Here, we focus on fragmentation of a primordial gas cloud collapsing in a 
CDM halo, with properties as found in our cosmological simulation. In our simulation, we find that the central gas cloud 
does not break up into multiple objects but remains intact as a single entity until the last output time. To study cloud 
fragmentation in detail and more generally, we examine the stability of the gas against thermal instability. 

A well-known stability criterion is the so-called Field criterion: 

where L is the energy loss rate per unit mass (Field 1965). Because the molecular hydrogen cooling rate is density- 
independent at high densities, and has a steep temperature dependence oc T" with a > 1 at 100 < T[K\ < 10000 (see Fig. 
1), the Field criterion is satisfied for primordial gas in this region for a fixed molecular fraction. Thermal instability can 
also be driven by a change in chemical composition, such as ionization/recombination (Defouew 1970; Goldsmith 1970) 
or molecule formation/dissociation (Yoneyama 1973). Sabano & Yoshii (1977) and Silk (1983) suggest that primordial 
gas can be chemo-thermally unstable. Chemo-thermal instability can be triggered when a rapid increase in the coolant 
fraction (H2 molecules) induces efficient cooling and condensation, leading to an enhanced production of molecules via 
three-body reactions. There are two proposed regimes where the chemo-thermal instability can be triggered: one at 
n ~ 10"'^°cm~^ where three-body reactions make the overall gas cooling very efficient, and the other at rt ~ lO^^cm"'^ 
where collision-induced emission cooling comes into play. 

Previously, ABN02 concluded that gas clouds do not fragment because of efficient turbulent mixing. Using one-zone 
calculations, Omukai & Yoshii (2003) conclude that the cloud core becomes unstable, but does not fragment because the 
cloud is already compact. Ripamonti & Abel (2004) study the thermal instability in the high density regime and argue 
that fragmentation is unlikely to occur. None of these previous works provide definitive answers, unfortunately. The 
high-resolution simulation of ABN02 employs the optically-thin approximation for H2 cooling, and thus their results do 
not provide accurate estimates for relevant quantities at densities above n > lO-'^^cm"'', whereas the other studies based on 
zero- or one-dimensional calculations suffer from the fact that fragmentation is intrinsically a three-dimensional problem. 
We are able to, for the first time, directly address this issue using a three-dimensional simulation with sufficient mass- 
and "physical" resolution in a fully cosmological context. 

We follow Omukai & Yoshii (2003) to investigate the stability of the gas cloud core against isobaric perturbations. 
For perturbations in temperature, density, and molecular fraction, ST,Sp,Sx cx exp(wi), under the constraint SP = 
5T + Sp + Sx = 0, linear stability analysis yields the dispersion relation 



Auj^ + Buj + C^ 0, 



(48) 
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B = ^{TLt - pL, -L)- ^{TFt -pF,-F)-A U + ^pF, + ^) , (50) 
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Fig. 7. — The instability growth parameter Q as a function of gas density. The shaded contour shows the distribution of the cloud core gas 
at a time when the central density is ^ lO^^cm""^. The thick solid line shows the evolution of the most dense particle in this Q — n plane. 
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In the above expressions, / is the molecular fraction (/ = 1 for a fully molecular gas), p ~ 2/(2 — /), F = df /dt is the 
net formation rate of hydrogen molecules, and L is the rate of energy loss per unit mass. Subscripts denote the respective 
partial derivatives. We have defined the local dynamical time 
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The dispersion equation H48|) has a positive real root (a growing mode solution) if and only if C < 0, because A is 
positive by definition. Furthermore, for perturbations to grow in a gravitationally collapsing gas, the characteristic growth 
timescale tg on 1/lu must be shorter than the dynamical time tdyn- Hence, we define the growth parameter 
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(53) 



If Q > 1, perturbations are expected to grow faster than the gravitational contraction of the cloud (Omukai & Yoshii 
2003; Ripamonti & Abel 2004). 

Fig. [7| shows the growth parameter evaluated locally for all the gas elements at the time when the central density is 
Tic = 10"'^^cm~^. We also show the evolutionary track for the cloud core by the solid line. We see that a part of the 
gas cloud is in the unstable region where Q > 1. However, Q never becomes much larger than unity. It is indeed less 
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than 1.5 for almost all the gas. Omukai & Yoshii (2003) argue that Q must be significantly larger than unity for a gas 
parcel to break up into multiple objects. This is because the size of the cloud, 1^ ~ Cgidym must be much larger than 
the length of growing perturbations l-p ~ Cstg. In other words, while a gas parcel can become thermally unstable, for the 
gas to break up into multiple objects, Q must be at least larger than two so that the perturbations can grow before the 
region gravitationally contracts. Fig. [7| shows that the cloud center never enters the regime Q > 2. Our results therefore 
support the conjecture of Omukai & Yoshii that the cloud becomes chemo-thermally unstable but evolves into a single 
object, because the cloud is already too compact. 

We do the same analysis at the time when the cloud core enters the high-density regime. Evaluating equations H48|l - H51|l . 
we find that the gas is always stable. To see this more clearly, we show the gas distribution in a thermodynamic phase 
diagram in Fig. |S1 We indicate the region where Q > by grey contours. For estimating the value of Q as a function of 
density and temperature, we calculate /hj assuming chemical equilibrium. We solve the cubic equation for the equilibrium 
molecular fraction 

fc22(l - ffr? + ifc23/(l - Ifrv' = hufil - fW + \k2ifn\ (54) 

where ki denote respective reaction rates tabulated in the Appendix. While the assumption of chemical equilibrium 
does not strictly hold in all the plotted range of density and temperature, it is indeed a good approximation for the 
relevant region where the actual gas evolutionary track lies. As clearly seen in Fig. |H1 the gas never approaches an 
unstable region. It appears that, since the cooling rate is larger at higher density and/or at higher temperatures (see 
the temperature dependence of CIE cooling in Fig. |2Jl, the gas evolutionary track circumvents the instability region. 
It is implausible that a gas parcel enters the instability region from the lower temperature side, because the cause of 
the instability is radiative cooling, which always tend to bring the gas temperature downward. We further argue that, 
because dissociation of hydrogen molecules occurring at T > 2000K effectively absorbs the heat input by contraction and 
thermalization, the gas evolves approximately isothermally around T ^ 2000K until full scale dissociation is completed, 
and the temperature cannot increase much without the density increasing. We therefore conclude that fragmentation 
owing to thermal instability does not take place in this regime either. 
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Fig. 8. — Stability analysis in the high-density regime. We plot the distribution of gas in the cloud core in the temperature-density phase 
space. The shaded regions indicate regions with Q = 0.1,0.5, 1 from light to dark grey, respectively. The gas does not enter the instability 
region. 



5.3. Stability against gravitational deformation 

Gas fragmentation could be triggered if a cloud is significantly flattened (Tohline 1981; Miyama et al. 1984; Tsuribe & 
Inutsuka 1999; Bromm, Coppi, Larson 1999) or if it contracts into filamentary structures (Uehara et al. 1996; Nakamura 
& Umemura 2001). Although Bromm et al. (1999) present a particular case which yields a disk-structure, their simulation 
is set up as such, by imposing an initial spin. In cosmological simulations, primordial gas clouds with such shapes are 
not generally found (Abel et al. 2002; Yoshida et al. 2003). There are two important mechanisms that lead to cloud 
fragmentation; one is the growth of deformation during collapse and the other is rotation-induced fragmentation. While 
these two mechanisms are, in general, coupled to each other, we here examine each effect separately. We first show that 
the dense, collapsing gas cloud formed with nn ~ lO^'cm"'^ would not fragment again at higher densities. The collapse 
of the dense core proceeds in a self-similar manner. This so-called Larson-Penston-type self-similar solution is known 
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to be unstable to non-spherical deformation if 7 = dlogP/dlogp is less than the critical value, 7crit = 1-097 (Hanawa 
& Matsumoto 2000; Lai 2000). An unstable core elongates by this instability and eventually fragments. Let us define 
the core's elongation as £ = {b — a) /a, where a and b are short and long axis lengths, respectively. The growth rate of 
the elongation v — d\og£/d\ogp is calculated by Hanawa & Matsumoto (2000) and Lai (2000) using Hnear perturbation 
theory. 
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Fig. 9. — Linear stability analysis for deformation growth. We plot the predicted evolution of the elongation ratio £ for the core (solid line). 
The evolutionary track of the core stays below S = 1 after the first run-away collapse at riu ~ lO'^cm"^, suggesting that the core is stable 
against deformation to filamentary structure. The dashed line shows the evolution of the effective ratio of specific heats, and the dotted line 
is the critical value 7crit = 1-097 from linear theory. 

Using the simulation outputs, we calculate the effective ratio of specific heats 7 = dlogPc/c'logPc at the density maximum 
as a function of density- The evolution of the core elongation is then calculated as 



£ = £q exp / '^('t-h) dlnn^ (55) 

for the initial elongation £0 at nn.o- In Fig- HI we show the evolution as a function of the central number density. The 
initial state is taken at the maximum of £ around lO^cm"'^, where the first run-away collapse begins. Early in the evolution 
(tt-h < lO^cm"'^), 7 exceeds 7crit, and the elongation decays; the core becomes rounder. With the onset of three-body H2 
formation, 7 starts decreasing and falls slightly below 7crit- There, the elongation grows but only slowly As is seen in 
Fig. 1^ it at most recovers the initial value £q. Since the effective 7 gradually increases again at n > 10^°cm~'^ when H2 
line cooling becomes inefficient, the deformation growth rate decreases. We thus conclude that the core must be stable 
initially and remains stable to this mode of fragmentation throughout the collapse, as is indeed the case in our simulation 
(see Fig. (S)). 

Next, we examine the effects of rotation. Tohline (1981) and many subsequent works (e.g. Tsuribe & Inutsuka 1999) 
conclude that basically two parameters, the cloud's initial ratio of thermal to gravitational energy (ap) and the initial 
ratio of rotational to gravitational energy (/3o), determine whether or not the cloud will eventually fragment, for a given 
adiabatic exponent 7. Here, the parameters are defined as, respectively, 

5c2i?o . ^IRI .... 
"°=2GM' ^°=3GM- 

Note that the latter quantity is essentially the degree of rotation support (see Section below). Tsuribe (2002) investi- 
gated a particular case for 7 = 1.1, which is close to the actual value for primordial gas, and concluded that the stability 
criterion is roughly given by > 0.3 for < /So < 0.3. (Note that the stability region does not differ much for a gas cloud 
with a constant initial density or for a gas cloud with central condensation [Tsuribe 2002].) We calculate these quantities 
at the time when the first run-away collapse begins, namely when the core density is ~ lO^cm^"^ and the temperature is 
^ 200K (Point C in Fig. E)). We find ac = 1.0 and (3c — 0.1, which satisfies the stability condition. In summary, the gas 
cloud in our cosmological simulation is found to be stable against deformation and rotation-induced fragmentation. We 
have validated the numerical results using analytical stability criteria. 
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Fig. 10. — Instantaneous gas mass accretion rate around the proto-stellar core (solid curve). We evaluate equation I57i and plot it as 
a function of enclosed gas mass. For reference, the dashed line shows the accretion rate obtained by Omukai &; Nishi (1998) using their 
post-singularity, self-similar solution. The dotted line is the critical mass accretion rate given by equation 1651 . 



5.4. Gas mass accretion rate 

We have seen that the prestehar cloud does not fragment into niuhiple clumps. Since the entire cloud has a much 
larger mass than the proto-stellar "seed" , the protostar will evolve in an inside-out manner, by accreting a large quantity 
of surrounding gas. The rate of accretion of the infalling gas is among the most important quantities in proto-stellar 
evolution. Since the late time evolution of the protostar itself affects the matter accretion rate, we are not able to measure 
the true mass accretion rate onto the protostar directly. Nevertheless, the instantaneous accretion rate likely provides a 
good estimate. In Fig. 1101 we show the instantaneous gas mass accretion rate at the last output time 

M ^4Trpr^v,^d{r) (57) 

using again the mass coordinate. We compare our result with that of ON98 which is derived from the post-singularity 
self-similar solution (Yahil 1983). The model accretion rate of ON98 has a shallower slope (i.e., larger accretion rate in the 
outer regions) than that of our simulation. This is reasonable because the ON98 model is based on the Larson-Penston 
type self-similar solution that predicts higher infall velocities at large distances because of the assumption of the solution 
itself; the velocity difference owes to the initial conditions and boundary effects, which are self-consistently calculated in 
our simulation. 

A simple, order of magnitude estimate for the mass accretion rate is given as a function of sound speed as 

A/-§, (58) 

and thus it is large for a high temperature gas. For T ^ 1000 K, the above estimate yields lO^'^AfQyr"^, similar to the 
simulation result. It is illustrative to compare this value with that in present-day star-forming regions which have T 10 
K. The temperature difference naturally explains the much larger accretion rate, by more than two orders of magnitude, for 
the primordial case. Intriguingly, the overall feature of the accretion rate shown in Fig. 1101 can be explained qualitatively 
from the temperature profile. The gradual increase at Monc = 100 — lGGGAf0 corresponds to the gradual temperature 
increase shown in the top-right panel in Fig. |B1 The slight shallowing at Monc ~ IOM0 corresponds to the temperature 
"dip" where three-body reactions promote molecule formation. Inside the radius, Mcnc < IOMq, the temperature rises 
again, and the accretion rate increases as well. The decline at the inner-most part is simply because of the output timing. 
At this time, the radial velocity gradually decreases towards the center (Fig. O; i.e. a radiative shock has not yet formed. 
In the outer part, Menc > IOOOMq, the increasing accretion rate owes to gravitational pull by the host dark matter halo, 
which is not seen in cloud evolution calculations without dark matter. 

The instantaneous accretion rate is larger than ABN02 found over the entire plotted mass range. While the small 
accretion rate at Monc < lOM© in ABN02 can be understood as a result of over-cooling (see Fig. 2 in ABNG2), it is unclear 
why the large-scale (A/onc ^ IOOMq) mass accretion rates differ. It cannot be attributed to the degree of rotational 
support because we obtained quite similar rotation values to ABN02 (see the next section). To study this discrepancy 
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further, we carry out two controlled simulations by embedding a primordial gas within an NFW dark matter halo similarly 
to the simulation presented in Sectional For one case we set /b — 0.16, consistent with our ACDM cosmology, whereas 
for the other case, we set /b = 0.05, that is the value adopted in the standard cold dark matter simulation with Jim = 1 
of ABN02. For both cases we assign zero velocities initially and let the gas cloud collapse in the halo's potential with the 
same mass. We measured the resulting mass accretion rates at the time when the central density reaches — lO^^cm"^. 
We find, interestingly, that the low baryon fraction case agrees well with the result of ABN02. The overall cooling time 
in the low baryon fraction case is longer than in the other one and hence collapse takes place more slowly. In the high 
baryon fraction case, external pressure from the ambient gas at the onset of collapse may also affect the accretion rate 
(Hennebelle et al. 2004; Motoyama & Yoshida 2003). Given the insignificant difference of up to a factor of a few, we do 
not discuss the discrepancy further in the present paper. Recently, O'Shea & Norman (2006a, b) performed a series of 
cosmological simulations to examine various environmental effects, such as large-scale local density and halo formation 
epoch, on the instantaneous gas mass accretion rate. They indeed find a substantial variation in accretion rate among 
different halos. We mention that very different cosmological environments, such as those realized in the early large-scale 
structure simulation of Gao et al. (2005) could have a noticeable impact on the mass accretion rate. 



The gas cloud starts collapsing with a finite amount of angular momentum (AM), and thus it is expected to spin- up 
gradually as it contracts, unless there are some mechanisms to transfer AM. Rotational support of the gas may play an 
important role in the process of accretion. To first order, we expect that mass accretion continues only if rotation does 
not halt collapse; i.e. if the angular momentum of the gas is small. We calculate the specific angular momentum profile 
around the protostar at the final output time. Fig. Illl shows the profile as a function of enclosed gas mass. The power-law 
profile is typical for gravitationally collapsed objects (Bullock et al. 2001; van dens Bosch et al. 2002). The AM profile's 
shape and even the amplitude are quite similar to the one found in the simulation of ABN02 over the entire plotted range. 
The degree of rotational support defined as 



is shown by the dashed line. Here, Lgp is the specific AM of the gas within radius r, and Vkcp is the Kepler velocity at 
that radius. It is about 40-50% for Menc = 0.01 — lOOOM©, and rotation is not (yet) able to stop the gas inflow. We 
have also measured /r at previous output times and found that /r has been monotonically increasing at small radii. To 
understand the reason for the characteristic AM profile, i.e., why low angular momentum material resides in the central 
regions, we examine two quantities. 

First, we make a conjecture that the region with little angular momentum collapses first. By measuring the probability 
distribution of specific AM in the cloud before the onset of collapse, we show that this is indeed the case in our simulation. 
At a time when the maximum density reaches lO^cm^^, we define the "core" as the region where the density is larger than 
half that of the center. The center is defined to be the position of the most dense particle. We then calculate the specific 
angular momentum with respect to random points within the core, by taking a fluid parcel of mass 0.1Mjcans,coic = 30Mq 
around the points. The resulting probability distribution is shown in Fig. 1121 There, the arrow indicates the value for the 
"true" center that collapses first and becomes the densest part in the later output time. We interpret this, together with 
the specific AM distribution, as evidence that regions with low angular momentum collapse fast in the (weakly-)turbulent 
medium. There is alway a central core within which density fluctuations are erased by soundwaves (hence having a size 

Cgt), and there is also an AM distribution within it. 

We further examine the evolution of AM in detail, making use of the Lagrangian nature of the SPH method. We trace 
the inner-most particles and calculate the evolution of the specific AM. We mark gas particles according to their mass 
coordinate at the final output time and divide them into four shells. We then keep track of the positions and velocities 
of the marked particles at earlier output times. At each output time, we shift the spatial coordinate such that the "true" 
central particle is always at the origin of the coordinate. (We note that the specific AM profile varies considerably if 
we define the center as the most dense region at each output time.) Fig. ^] shows the evolution of the specific angular 
momentum for the four mass shells. As is clearly seen in the right panel of Fig. 1121 the specific angular momentum of 
each mass shell is well-conserved, while the mean radii decreased by a factor of 3-5. This causes the gradual increase of 
rotational support in the inner region as discussed previously. 



5.5. Angular momentum 




(59) 
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Fig. 1 1. — Specific angular momentum profile in the gas cloud. The dashed line shows the degree of rotational support as defined by 
equation 1591 . It is about 40-50% for the entire cloud. 

In order to verify that numerical viscosity did not significantly influence our results, we have explicitly checked the 
torque exerted on the gas cloud. In SPH, the equation of motion for each gas particle is expressed as 



dv 

m- = F,,., 



pressure 



(60) 



where the last term denotes the force owing to artificial viscosity (see Springel 2005 for exact expressions of each term). 
Then the torque N = dL/dt exerted on a gas cloud consisting of M particles can be reduced to a sum of three components: 



M 
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We store the forces on each particle in output files and calculate explicitly all the terms in equation H61I) . We find that the 



dominant torques are A'grav and 
O.lpc), A^grav and iV, 



^ pressure 



throughout the evolution. Within the collapsing gas cloud (Mone < 2OOM0, 



* pressure 



are completely dominant, and Advise is always less than 10% of iV, 



pressure 



r < 

for all the mass shells. 

Note that a small amount of dissipation by weak shocks is expected, because there are some fiuid elements that are 
moving at transonic velocities, as found in the adaptive mesh refinement simulations of ABN02. We have also performed 
the spherical collapse test of Norman, Wilson & Barton (1980). We set up a rotating spherical cloud as in Norman et 
al. (see also Truelove et al. 1998); M = 1Mq,R = 7.01 x lO^^cm and the initial rotation n ^ 3.04 x lO^^^rad sec'^ 
We use a half million particles to represent the sphere. For an ideal (inviscid) fluid with no means of angular momentum 
redistribution, the specific AM distribution is conserved, and it takes a simple analytic form 



M(< K) = M 



K 



3/2' 



(62) 



where M(< K) is the mass of fluid elements that have speciflc AM less than K. We measure the distribution at the 
initial time and that after about one free-fall time. The AM distribution is found to be conserved rather well, with the 
deviation at the smallest mass scale being less than 20%. This is a similar level of agreement as found in, for example, the 
calculation of Truelove et al. (1998). Although a good level of AM conservation in this test is expected for a Lagrangian 
hydrodynamics code, our test results reassure that the artificial viscosity does not do any harm to angular momentum 
transport in dynamically collapsing gas spheres such as those studied in the present paper. 
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Fig. 12. — (Left) Probability distribution of specific angular momentum in the core region when the first collapse is being triggered. The 
part that becomes the central region in the end has a small angular momentum initially (indicated by the solid arrow), whereas the temporal 
center (most dense region) before the collapse has a larger angular momentum of ~ 0.03 pc km/sec. (Right) Evolution of the 'Lagrangian' 
specific angular momentum with respect to the marked center of the gas cloud for four groups of gas elements that are labeled according to 
their final mass coordinates. 



5.6. Evolution of the protostar 

It is still beyond the capability of the current simulation code and architecture to directly follow the evolution of the 
protostar to the point where it settles onto the zero-age main sequence (ZAMS). We tackle this demanding problem by 
employing the scheme for protostellar evolution of Stabler, Shu, & Taam (1980a, b) as modified by Stabler et al. (1986a,b) 
and Omukai & Palla (2001; 2003). We assume that the gas accretion takes place in an approximately spherically symmetric 
manner. We use the mass accretion rate shown in Fig. ^| as an input to the protostellar evolution code. While this 
assumption may seem over-simplified, the gas mass accretion rate from our cosmological simulation provides a reasonable 
estimate for true accretion rates onto a protostar, at least for the initial accretion phase, because the central gas cloud 
is roughly spherical and the gas in the envelope is not yet in a disk structure. Moreover, the protostellar evolution 
calculations include various physics and thus will give a more self-consistent solution of the evolution of a protostar for a 
given mass accretion rate. We thus employ this method rather than use simple time-scale arguments as done by ABN02 
and Bromm & Loeb (2004). 

Briefly, the evolution of a protostar is treated as a sequence of a growing hydrostatic core with an accreting envelope. 
The core is assumed to be in hydrostatic equilibrium, and the ordinary stellar structure equations are applied. The 
structure of the accreting envelope is calculated with the assumption that the flow is steady for a given mass accretion 
rate. Only the inner envelope that is optically thick to continuum is included. The effects of radiation pressure are 
considered in constructing the structure of the envelope. The details of the calculation are found in Stabler, Palla, & 
Salpeter (1986a,b) and Omukai & Palla (2001; 2003). 

We use the time-(equivalently, protostellar mass-) dependent accretion rate Mfid derived from our simulation (Fig. I1U|) . 
We refer to this rate as the fiducial rate hereafter. The evolution of the protostellar radius is shown in Fig. E| as a function 
of the protostellar mass. Considering the possibility that some fraction of infalling matter stacks in the circumstellar disk 
or is ejected as an outflow, we also show cases of accretion rates reduced from our fiducial value by a factor of 2/3 and 

Initially, the Kelvin-Helmholz (KH) time t^ii, which is the timescale for a star to lose its thermal energy by radiation, is 
longer than the accretion timescale tacc because the temperature of the protostar is low and its opacity owing to free-free 
absorption is large {ne oc T"^/^). The accreted material simply piles up on the stellar surface without cooling. This phase 
is called the adiabatic accretion phase and lasts up to A/* ~ IOMq. During this phase, the radius first increases linearly 
(in log scale) with the protostellar mass, and then remains almost constant at 100i?Q. This behavior can be explained as 
follows. 

The protostellar radius in the adiabatic accretion phase depends on the protostellar mass and accretion rate as (Stabler 
et al. 1986a,b) 

oc M^'^M^-^. (63) 

The increase of the radius in the earliest phase of evolution is explained by the increasing mass and accretion rate shown 
in Fig. 1101 This is partly because our mass resolution is limited, and hence the output time is slightly before the formation 
of an accretion shock. Also, adjustment from a somewhat arbitrary initial structure to the appropriate accretion phase 
produces somewhat complex behavior at il/* < O.3M0. The accretion rate decays with mass as cx M^'^'^ in the range 
IMq < < IOMq in our fiducial model. Substituting this into the above relation, we see that the mass dependence of 
the radius almost cancels. This explains the approximate constancy of the radius in this range. 
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Fig. 13. — Proto-stellar evolution. The assumed accretion rates are taken from our fiducial model Afflj (solid), (2/3)Mfl(j (dashed), 
(l/3)Mfi;i (dash-dotted). The solid points indicate the time when hydrogen burning begins. 

With increasing protostellar mass, increasing temperature in the stellar interior makes the opacity lower and the radi- 
ation diffuses out more easily. When the protostellar mass reaches ~ IOMq, heat deposited in the interior is transported 
outwards by a luminosity wave, whose arrival to the surface makes the star swell suddenly (Stahler et al. f 986a,b). After 
that, the protostar enters the Kelvin-Helmholtz contraction phase. By losing thermal energy to radiation, the protostar 
shrinks as (Omukai & Palla 2003) 

cx M.MZ'^. (64) 

The interior temperature increases until the contraction is halted by nuclear burning. 

Energy generation by the p-p chain is not sufficient to stop the contraction of the already massive star. When the 
central temperature reaches lO^K, hydrogen burning by the CNO cycle begins with a slight amount of C synthesized 
by He burning. The onset of hydrogen burning by the CNO cycle is marked by a solid circle in the figure. The energy 
generation by hydrogen burning halts contraction around IOOMq, and the star reaches the ZAMS. The protostar relaxes 
to a ZAMS star eventually within about 10^ years from the birth of the protostellar seed. In our calculation, there is no 
sign that accretion is halted by stellar feedback although the radiation pressure onto the accreting matter is included. 
Even after its arrival onto the ZAMS, the star continues to accrete. In principle, the star can grow in mass during its 
entire lifetime. (Note that the accretion rate at the outer boundary is given as a boundary condition. The gas accretion is 
sustained because of this condition. See the discussion below.) The stellar mass would become greater than a few hundred 
solar masses during its main sequence lifetime (~ 3 x 10^ yrs) if we were to allow mass accretion to continue using our 
estimated accretion rate. This mass scale is similar to the total gas mass of the dense core. This coincidence is not by 
chance: since the stellar lifetime is longer than the free-fall time of the dense core (~ 3 x 10^ yrs for nn ^ lO^cm"'^), the 
whole core falls onto the star. 

In order for the accretion to continue, the total luminosity of the protostar Ltot, given by the sum of the luminosity 
from the stellar interior and that from the accretion shock Lace — GM^M^/ R^, must not exceed the Eddington 
limit iEdd — 47rcG'M»/Kos, where k^s is the electron scattering opacity, the dominant source of opacity in the inner 
envelope. The total luminosity approaches the Eddington limit as the star approaches the ZAMS, decreasing the radius 
and increasing the luminosity. To reach the ZAMS, the mass accretion rate must be less than 
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(65) 



where i?zAMS a-nd Lzams £^re the radius and luminosity of the ZAMS star of mass M^. The dependence of the critical 
accretion rate on stellar mass is very weak and becomes approximately constant at Merit — 4 x lO~'^M0yr~^ (Omukai & 
Palla 2003). This means that once the star reaches the ZAMS, the luminosity remains always below the Eddington limit 
and the accretion continues unimpeded. In Fig. 1101 we show the critical rate by the dotted line. Although the accretion 
rate is higher than the critical value in the early evolution, the total luminosity falls below the Eddington limit owing 
to the larger radius and smaller interior luminosity in that phase than a ZAMS star of the same mass. For the fiducial 
accretion rate, the accretion rate drops below the critical value at M* ~ 3OA/0, well before the star reaches the ZAMS. 
Consequently, the total luminosity never exceeds the Eddington limit and accretion continues. 
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Wc also explored two other cases besides our fiducial model. For these cases, we attempt to model crudely the feedback 
effects from the protostar and the consequences of rotation by reducing the mass accretion rate from the original value. 
The formed protostar could launch a protostellar outflow if a small magnetic field is present initially (Machida et al. 
2006). Some fraction of infalling matter in the envelope is eventually ejected as an outflow without falling onto the central 
star and the accretion rate onto the protostar is reduced. To this end, we use reduced accretion rates of (2/3)Mfid and 
(l/3)Mfid- The overall evolution of these cases is similar to that of the fiducial case, but characteristic mass scales are 
shifted to somewhat lower values. For the case with M = (l/3)Mfid, the protostar reaches the ZAMS when its mass is 
~ 6QMq. In the case of disk accretion, not only the accretion rate, but also the surface boimdary condition of the star is 
altered. For smaller accretion rates (10~^ — 10~^M©/yr), Palla & Stabler (1992) studied the cases of different boundary 
conditions, which correspond to spherical or disk accretion cases. Intriguingly, they showed that the evolutionary behavior 
is quite similar in both cases qualitatively; the star reaches ZAMS somewhat earlier in the disk accretion case. 

In summary, the accretion is not halted by radiative feedback for the fiducial and reduced accretion rates. Note, however, 
that wc only considered the inner envelope, whic;li is optically thick to continuum. In the outer envelope, the Lya opacity 
could be important for halting accretion (Doroshkevich & Kolesnik 1976; Tan & McKee 2004). The expansion of an Hii 
region might be an another possibility to hinder accretion (Larson & Starrfield 1971; Tan & McKee 2004). Since the 
expansion of an Hil region is quenched for a large accretion rate assuming spherical symmetry (e.g., Omiikai & Inutsuka 
2002), non-spherical effects must be included for evaluating the exact time behavior of the expansion of the Hii region. 
While these uncertainties remain in determining the exact mass-scale of the first stars, the fact that they are rather 
massive seems to be robust. 

6. SUMMARY AND DISCUSSION 

We have developed numerical techniques to implement the required chemical and radiative processes for following 
the thermal evolution of primordial gases to very high densities. Using a cosmological simulation with hydrodynamics 
and chemistry, we have studied the formation and evolution of pre-stellar gas clouds in a ACDM universe. Our three- 
dimensional simulation for the first time resolves the inner ~ lAf© fully molecular core, allowing us to obtain correct 
density, temperature, and velocity structure around the primordial protostar. These quantities altogether yield an accurate 
gas mass accretion rate around the protostar that was uncertain in previous works. We have derived three important 
facts from our simulations: (1) the primordial gas cloud in a cosmological minihalo does not fragment, but yields a single 
proto-stcllar seed, (2) the cloud core has a small angular momentum so that rotation does not halt the collapse, nor is a 
disk formed, and (3) the rate of accretion from the cloud envelope is large. From these facts, we conclude that the first 
stars are massive. By performing a proto-stellar evolution calculation, we derive that the protostar grows to ~ lOOM© for 
the particular case found in our simulation. We derive the stellar mass for the first time from a self-consistent calculation 
of proto-stellar evolution, coupled with a cosmological simulation. 

While wc simulated only a single halo in detail, the particular case shows a convincing and plausible case for the 
formation of massive stars in early mini-halos in the CDM model. The fate of such massive primordial stars (Population 
III), whether or not they are commonly born in various environments, arc of considerable interest. They are likely 
responsible for various feedback effects in the early Universe (Haiman, Rees, & Loeb 1997; Bromm, Yoshida, & Hernquist 
2003; Yoshida, Bromm, & Hernquist 2004; see Ciardi & Ferrara 2005 for a review). Stars with mass 140 — 26OM0 are 
believed to trigger pair-instability supernovae and expel all the processed heavy elements in the explosion (Bond, Arnett, 
& Carr 1984; Heger & Woosley 2002). It has also been suggested that less massive stars die as hypernovae (Umeda &: 
Nomoto 2002), producing a quite different metal yield from the pair- instability case. Interestingly, the observed abundance 
pattern of extremely metal poor stars appears to favor the latter scenario (Umeda & Nomoto 2003; Frebcl et al. 2005). 
Energetic supernovae are also destructive, being effective in evacuating the halo gas (Bromm, Yoshida & Hernquist 2003; 
Wada & Venkatesan 2003; Kitayama & Yoshida 2005) and possibly quenching further star-formation in that region. 
A single supernova would suffice to pollute a large volume of the interstellar matter from which very metal-poor stars 
might have been formed. Metal-enrichment by the first supernovae may also lead to the formation of ordinary stellar 
populations including low-mass stars (e.g. Mackcy et al. 2003) if metal-enrichment is relatively confined within a small 
region. Intriguingly, Jimenez & Haiman (2006) argue that metal-mixing could be inefficient and that primordial stars are 
formed in galaxies at 2 < 5. 

Massive stars emit a large number of ionizing photons and can ionize a large volume of the intergalactic medium 
around them (Kitayama et al. 2004; Whalen et al. 2004). Recently, Yoshida (2006) carried out radiation-hydrodynamic 

calculations of early Hii regions and concluded that radiation from a massive Population HI star can completely ionize the 
gas within the host halo, quenching further gas condensation and star-formation for a significant period of cosmic time. 
Both massive and very massive (> 3OOA'/0) stars likely leave a black hole remnant. It is plausible that the remnants become 
the seeds for supermassive black holes that power luminous quasars found at z ^ 6 (Li et al. 2006, in preparation) and 
subsequently at lower rcdshifts (e.g. Di Matteo et al. 2005; Hopkins et al. 2005, 2006). Because bottom-up hierarchical 
structure formation is a generic prediction of the CDM model, it is natural to think that the early population stars affect 
and even set the initial conditions for primeval galaxy formation. It is important to explore feedback from the first stars 
in more detail under a proper cosmological set-up. 

The prospects for detecting the first generation of stars appear promising. Supernovae resulting from the collapse of 
massive stars can be observed if they trigger gamma-ray bursts (GRBs). There is growing evidence that long duration 
GRBs are associated with energetic supernovae (e.g. Hjorth et al. 2003). Gou et al. (2004) suggest that the afterglow 
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of high-redshift GRBs can be detected in X-rays. There are also a broad range of observations through which we can, 
in principle, indirectly observe star-formation at high redshifts, through, for example, angular fluctuations in redshifted 
21 cm emission from intergalactic gas as stars reionize the Universe (see, e.g. Zaldarriaga et al. 2004; Furlanetto et al. 
2004; Zahn et al. 2006), or from secondary anisotropics imprinted from stellar radiation on the CMB (e.g. McQuinn et 
al. 2005; Zahn et al. 2005). 

Finally, we remark that it is still too early to definitely conclude the exact mass or mass range of the first stars. The 
actual mass accretion could take place in a more complicated manner than we assume. Feedback from the protostar itself 
can affect the dynamics and thermal evolution of the accreting gas. Ultra-violet radiation from the protostar likely affects 
the gas infall rate. However, primordial gas has a much smaller opacity because of the absence of dust, and thus radiation 
pressure is generally not important (Omukai & Inutsuka 2002). The effect of radiation pressure can also be reduced if 
protostellar outflows provide optically thin cavities through which photons can escape (Krumholz, McKee & Klein 2005). 
Future calculations of proto-stellar evolution need to include these feedback effects onto the infalling gas. Also, accurate 
modeling of disk accretion may be necessary in general cases (Tan & McKee 2004). Although our simulation is currently 
resolution-limited (rather than physics-limited), with the maximum density being still five orders of magnitude smaller 
than stellar densities, we foresee that it will be possible to carry out an ab initio calculation of the formation of primordial 
protostars in the near future. 
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APPENDIX 



Table 1 

REACTION RATE COEFFICIENTS 



Reactions 



Rate Coefficients (cm^s ^) 



Reference 



(1) 



(2) 



(3) 



(4) 
(5) 



(6) 
(7) 
(8) 
(9) 
(10) 

(11) 
(12) 



H + e ^ H+ + 2e 



H+ + , 



He + e ^ He+ + 2e 



He+ + e^lle + hu 



He+ + e ^ He++ + 2e 



He++ + e^Re+ + hiy 
H + e H" + /ii/ 
H- + H ^ H2 + e 
H + H+ ^ H J + /iz/ 



H+ + H ^ + H+ 



Ha + H ^ 3H 

H2 + H+ ^ H+ + H 



fci = exp[-32.71396786+ 13.536556 (InTe) 

-5.73932875 (Inrc)^ + 1.56315498 (InTe)^ 
-0.2877056 (InTo)^ + 0.0348255977 (InTe)^ 
-0.00263197617 (InTo)^ + 0.000111954395 (InTe)'' 
-2.03914985 X 10"*^' (InTc)^] 

fca = exp[-28.6130338- 0.72411256 InTe 

-0.0202604473 (InTe)^ - 0.00238086188 (InTe)^ 

-0.000321260521 (InTc)^ 

-0.0000142150291 (InTe)-"^ 

+4.98910892 x 10"^ (InTc)** 

+5.75561414 X 10"^ (InTc)^ 

-1.85676704 X 10"^ (InTe)^ 

-3.07113524 X lO"'^ (InTe)^] 

fcg = exp[-44.09864886+ 23.91596563 InTe 

-10.7532302 (InTo)^ + 3.05803875 (lnTo)=' 
-0.56851189 (InTe)^ + 0.0679539123 (Inre)^ 
-0.00500905610 (InTc)*^' + 0.000206723616 (InTe)^ 
-3.64916141 X 10"*^ (InTg)^] 

kir = 3.925 X io-i3t^-o-6353 

k4d = 1.544 X 10-9T-^-5exp(-48.596/Te) x [0.3 + exp(8.10/Te)] 

fcs = exp[-68.71040990+ 43.93347633 In To 

-18.4806699 (InTe)^ + 4.70162649 (Inre)^ 
-0.76924663 (InTe)'' + 0.08113042 (InTe)^ 
-0.00532402063 (InTe)*^ 
+0.000197570531 (InTe)^ 
-3.16558106 X 10"*^ (InTg)^) 

fee = 2 X A;2(Te/4) 

kr = 1.4 X IQ-^^T"-^^^ exp{-T/ 16200) 
ks = 4.0 X 

fcg = dex[-19.38- 1.523 logr+ 1.118(logT)2 - 0.1269(logT)3] 

kio = 6.0 X 10-10 
fit by reference 6 

fci2 = exp(-21237.15/r) x [-3.3232183 x 10'^ 

+3.3735382 x 10"^ (InT - 1.4491368 x 10""^ (InT)^ 
+3.4172805 X 10"^ (InT)^ - 4.7813720 x 10"^ (InT)* 
+3.9731542 x lO^^" (InT)^ - 1.8171411 x 10"" (InT)^ 
+3.5311932 X 10-13 (InT)^] 



1 
1 
2 
2 

2,3,4 

5 

6 
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Table 2 

REACTION RATE COEFFICIENTS (continued) 



Reactions Rate CoefBcients (cm'^s ^) Reference 

(13) H2+e^2H + e fcig = 3.73 x IQ-" ^0.1121 ^xp(^gg430/2^) 9 

(14) He+ + H ^ He + H+ + hi^ ku = 1.20 x lQ-^^{T/2,mf-'^^ 5, 10 

(15) He + H+ ^ He+ + H fcis = 1.26 x iq-^t-^-^^ exp(-127500/T) (T < 10^ K) 

fci5 = 4 X 10-37^4-74 \t > 10^ K) 5, 11 

(16) H-+e^H + 2e Aiie = exp[- 18.01849334+ 2.3608522 In Te 

-0.28274430 (InTe)^ + 0.0162331664 (Inre)^ 

-0.0336501203 (InTe)"* + 0.0117832978 (InTe)^ 

-0.00165619470 (Inrc)^ + 0.000106827520 (InTe)'' 

-2.63128581 x 10"^ (InTo)^] 1 

(17) H- + H -> 2H + e kn = exp[-20.37260896+ 1.13944933 InTe 

-0.14210135 (lnTe)2 + 0.0084644554 (InTe)^ 
-0.0014327641 (InTc)* + 0.00020122503 (Inre)^ 
+0.000086639632 (InTe)^ 
-0.000025850097 (InTef 
+2.4555012 X lO"*^ (In To)* 

-8.0683825 X 10-* (In Te)9] 1 



(18) 


H- +H+ 


^ 2H 


kis 


= 6.3 X 10"* + 5.7 X io-ex-o-5 _ 9 2 x lO-^^T^-^ + 4.4 x IQ-^^T 


2, 12 


(19) 


H- +H+ 


^H+ + e 


kig 
ki9 


= 4.0 X 10-4T-i-4exp(-15100.0/T) (T > 10^ K) 
= 1.0 X 10-8T-0-4 (T < 10" K) 


13 


(20) 


H+ + e^ 


► 2H 


k20 


= 2.0 X lO-^T-O-5 


2 


(21) 


H++H- 


^H + H2 


k2i 


= 5.0 X 10-^(T/100)-°-5 


1 


(22) 


3H^H2 


+ H 


k22 


= 5 X lO-^^T-i 


8 


(23) 


2H + H2 - 


2H2 


k2Z 


= fc22/8 


8 


(24) 


H2 + H2 - 


2H + H2 


k24 


= 8.125 X 10-*T-i/2cxp(-52000.0/T) 





x(1.0-exp(-6000.0/T)) 



Note.— We denote To for T in eV, otherwise T is in K. References are (1): Abel et al. (1997); (2): Galli & Palla 
(1998); (3): Ramaker & Peek (1976); (4): Stancil ct al. (1993); (5): Glover & Brand (2003); (6): Martin et al. (1996); 
(7): Savin et al. (2004a,b); (8): Palla, Stabler, & Salpeter (1983); (9): Stibbe & Tennyson (1999); (10): Zygelman et al. 
(1989); (11): Kimura et al. (1993); (12): Peterson et al. (1971); (13): Shapiro & Kang (1987) 



